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Abstract — The term ’’iterative method” refers to a wide 
range of techniques which use successive approximations to 
obtain more accurate solutions .In this paper an attempt to 
solve systems of linear equations of the form AX=b, where A is a 
known square and positive definite matrix. We dealt with two 
iterative methods namely stationary (Jacobi, Gauss-Seidel, SOR) 
and non-stationary (Conjugate Gradient, Preconditioned 
Conjugate Gradient).To achieve the required solutions more 
quickly we shall demonstrate algorithms for each of these 
methods .Then using Matlab language these algorithms are 
transformed and then used as iterative methods for solving these 
linear systems of linear equations. Moreover we compare the 
results and outputs of the various methods of solutions of a 
numerical example. The result of this paper the using of 
non-stationary methods is more accurate than stationary 
methods. These methods are recommended for similar situations 
which are arise in many important settings such as finite 
differences, finite element methods for solving partial 
differential equations and structural and circuit analysis . 


Index Terms — Matlab language, iterative method, Jacobi. 


I. Stationary iterative methods 

The stationary methods we deal with are Jacobi iteration 
method, Gauss-Seidel iteration method and SOR iteration 
method. 

A. Jacobi iteration method 

The Jacobi method is a method in linear algebra for 
determining the solutions of square systems of linear 
equations. It is one of the stationary iterative methods where 
the number of iterations is equal to the number of variables. 
Usually the Jacobi method is based on solving for every 
variable Xj of the vector of variables X T =(xi,x 2 ,....,x n ) , 
locally with respect to the other variables .One iteration of the 
method corresponds to solve every variable once .The 
resulting method is easy to understand and implement, but the 
convergence with respect to the iteration parameter k is slow. 

B. Description of Jacobi ’s method: - 

Consider a square system of n linear equations in n variables 
AX=b (1.1) 

where the coefficients matrix known A is 
A=(fl£ p j) Lj = 1(1) ?i 

The column matrix of unknown variables to be determined X 
isX t = (x v x 2 xf) 

and the column matrix of known constants b is 


b t = (fi L , V--- bn) 

The system of linear equations can be rewritten 
(D+R)X=b (1.2) 


where A=D+R , D =( 0 ^) i = l(l)n 
is the diagonal matrix D of A and R = L T U , 
where L and IT are strictly lower and 
strictly Upper matrix of A 

Therefore, if the inverse 0 -1 exists and Eqn.(1.2) can be 
written as 


X=B -1 (b-RX) (1.3) 

The Jacobi method is an iterative technique based on solving 
the left hand side of this expression for X using a previous 
value for X on the right hand side .Hence , Eqn.(1.3) can be 
rewritten in the following iterative form after k iterations as 
(b-R^C*-i)) ,k=l,2,--- (1.4) 

Rewriting Eqn.(1.4) in matrix form . We get by equating 
corresponding entries on both sides 


Cfc] if, fJt-0 \ 

f£ y 

V j*L / 


i=l,...,n ,k=l,2,.... (1.5) 

We observe that Eqn. (1.4) can be rewritten in the form 

X ik'i = TX (k-Vi + c k=i j2j , (1.6) 

where T=D _1 (-L-U) and C = D -1 b 

or equivalently as in the matrix form of the Jacobi iterative 

method 


X^=D -1 (— L— iflX^-^+D'b k=l,2,.... 
where T=D _1 (-L-U) and C = D _1 b 

In general the stopping criterion of an iterative method is 
to iterate until 


11x0*11 


<E 


for some prescribed tolerance E> 0 .For this purpose, any 
convenient norm can be used and usually the L,, norm , i.e. 


IIkmIL 


l <e,x- 


0. 


This means that if X :k 1 - is an approximation to X — ,then 
the absolute error is| |X 1 ^’ — X ,:k-1J | , and the relative error is 


ll* k: 'IL 


provided that 


X M 


O.This implies that a 


vector X" can approximate X to t significant digits if t is the 

1 1 X S " 1 1 

largest non negative integer for which — — - <5*10 _L 


C. The Matlab Program for JACOBI 

The Matlab Program for Jacobi’s Method with it’s 
Command Window is shown in Fig. (1.1) 
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D. Gauss-Seidel iteration method 

The Gauss-Seidel method is like the Jacobi method, 
except that it uses updated values as soon as they are available 
.In general, if the Jacobi method converges, the Gauss-Seidel 
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method will converge faster than the Jacobi method, though 
still relatively slowly. 

Now if We consider again the system in Eqn.(l.l) 

AX=b (1.1) 

and if We decompose A into a lower triangular component 
L* and a strictly upper triangular component U. Moreover if 
D is the diagonal component of A and L is the strictly lower 
component of A then 

A= L*+U 

Where 


L* — L+D 


Therefore the given system of linear equations in Eqn.(l.l) 
can be rewritten as: 


(L* +U)X=b 


Using this we get 

L*X=b-UX (1.7) 

The Gauss-Seidel method is an iterative technique that solves 
the left hand side of this expression for X, using previous 
values for X on the right hand side. Using Eqn.(1.7)iteratively 
,it can be written as: 

X^= L: 1 (b-UX' h -^),k=l,2,.... , (1.8) 

provided that Li 1 exists. Substituting for L* from above we 
get 

= (D 4- L) 1 (b)-(D+L) 1 UV h - V 


k=l,2,...., (1.9) 

where T = -((D + L) ')U and OCD + LJ^b 
Therefore Gauss-Seidel technique has the form 

X®= T X (k - 1 >+C,k=l,2,.... (1.10) 

However ,by taking advantage of the triangular forms 
of D,L,U and as in Jacobi's 

method ,the elements of X :k - can be computed sequentially 
by forward substitution using the following form of Eqn. 
(1.5) above: 


tu_i A - \ 

^ ^ V — Kj (k_D ) 


( 1 . 11 ) 


i=l,2,...,n , k=l,2,3, 

rib 

The computation of Xj' ■ uses only the previous elements of 
X :k_1 ‘ that have already been computed in advanced to 
iteration k.This means that , unlike the Jacobi method, only 
one storage vector is required as elements can be over 
rewritten as they are computed , which can considered as an 
advantageous for large problems .The computation for each 
element cannot be done in parallel . Furthermore, the values at 
each iteration are dependent on the order of the original 
equations. 


E. The Matlab Program for Gauss-Seidel iteration 
method 

The Matlab Program for Gauss-Seidel Method with it’s 
Command Window is shown in the Fig. (1.2) 


F. SOR iteration method (Successive Over Relaxation) 

SOR method is devised by applying an extrapolation w to 
the Gauss-Seidel method .This extrapolation takes the form of 
a weighted average between the previous iteration and the 
computed Gauss-Seidel iteration successively. 

If w>0 is a constant, the system of linear equations in 
Eqn. (1.1) can be written as 

(D+wL)X= wb- [ wU+( w- 1 )D] X (1.12) 

Using this Eqn.(1.12) the SOR iteration method is given by 


X Ck? = (D + WL) -1 [(1 - w)D — wU]X Ek - D 
+w(D -I- wO“ L b k=l,2,3... (1.13) 

provided 1 -. 3 -|- u-’Lj -1 exists . 


function j acobt2( A. b. tol. x. K) 

%jacobi(A. b. tol.X.N) solve iteratively a system of linear equations %whereby 
%A is tie coefficient matrix, and bis the right-hand side column vector 
%tol : relative residual error tolerance for break condition 
%X : Nxl start vector (the initial guess) 

%N is the maximum number of iterations 

%The method implemented is the Jacobi iterative 

%The starting vector is the null vector, but can be adjusted to one's needs 

%The iterative form is based on the Jacobi transition iteration matrix 


%Tj = mv(D)*(L+U) and the constant vector cj = mv(D)*b 


%The output is the solution vector X 


%splitting matrix A into the three matrices L U and D 
D = diag(diaig(A)); 

L = tril(-A-l); 

TJ= trm(-Al); 

%transition matrix and constant vector used for iterations 
Tj = inv(D)*(L+U); 


cj = inv(D)*b; 
k= 1; 

while k- r -=N 
x(:_k+l) = Tj * x(:,k) + cj ; 

B= norm(x(: ; k+l)-x(:,k)); 
if B^tol 

dispCIhe procedure was successful 1 ) 


» A=[Q.2 0.11 10; 0.1 4-11 -1; 
1-1600-2; 1 10 8 4; 
0-1-2 4 700] ; 
»x=[0;0;0;0;0]; 

» tol=0.00005 ; 

» N=91 ; 

»jacobi2(Ab,toLx,N) 

The procedure was successful 
Condition xTk+1) -x'\k) <tol was 
met after k iterations 
91 
x= 

7.8597 
0.4229 
-0.0736 
-0 5406 
00106 


dis pCC ondition ||x'Xk+l.) - x'Nk)|| * = tol was met after k iterations J 


dtspvk); dispCx = ") ;dis pf x( yk+ 1 )); 
break 


k = k+l; 
end 

if B-tol || k > N 

di 3 pv'U'.fa-v inmvrn numb er of iterations reached without satisfying co nditio n : "} ;di3 p£ w'Xk+l) - x"ik)| 
^tol r >; disputed); 

dispCPleaae, examine the sequence of iterates 1 ) 

dispfThcase you observe convergence, then increase the maximum number of iterations') 

dispCx); 

end 

Tig.(l.l) The JACOBI Matlab Pro gram with its 
Command Window for RUNNING it. 


function gauss_sddd(A, b,toi,x, N) 

%Gauss_Seidel(A b,toLX : N) solve iterativdy a system of linear %equationswhereby 
%A is the coefficient matrix, and b is the right-hand side column 
%vector 

%tol : relative residual error tolerance for break condition 
%X : Nx 1 start vector (the initial guess) 

%N is the maximum number ofiteratians 
%The method implemented is the Gauss-Seidel iterative 
%The starting vector is the null vector, but can be adjusted to one's needs 
%The iterative form is based on the Gauss-Seidd transition, iteration matrix 
%Tg = inv(D-L)*U and the constant vector eg = inv(D-L)*b 
%The output is the solution vector x 
%splitting matrix A into the three matrices L, U and D 


D = diag(diag(A)): 
L = tril(-A,-1); 

U = triu(-A,l); 


» A= [0.2 0.1 1 1 0 ; 0.1 4 -1 1-1:1-160 0-2; 

1 10 8 4:0-1-24 700]; 

» x=[0;0;0;0;0]; 

» b= [1;2;3;4;5]; 

» N=32; 

» tol=0.00005 ; 

» gauss_seidel(A,b,td,x,N) 

The procedure was successful 
Condition ||xXk + l) - x A (k)|| <tol 
was met after k iterations 3 1 


7.8596 

0.4229 

-0.0736 

-0.5406 

0.0106 


° otransition matrix and constant vector used for iterations 
Tg = inv(D-L)*U; 
cg = inv(D-L)*b: 
k= 1; 

while k<=N 
x(:)c+l) = Tg*x(:,k) + eg; 

B=norm(x(:,k+l)-x(:,k)); 

If B<tol 

disp('The procedure was successful’) 
dispCCondition Hx'Xk+l) - x-Xk)|| <tol 
was met after k iterations) 
disp(k); disp(’x = ');disp(x(:Jc+l)); 
break 
end 

k=k+l; 

end 

if B>tol || k>N 

disp('\laximum number of iterations readied without satisfying condition:) 
dispdlx'Xk+l) - x'Xk)|| <tol'); disp(td); 
disp('Please, examine the sequence of iterates') 

disp('In case you obsewe convergence, then increase the maximum number of iterations') 

disp(x); 

end 

% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&&&& 


% A=[0.2 0.1 1 1 0 : 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 : 0 -1 -2 4 700] 
% N=20 
% b=[l:2:3:4:5] 

% x=[0;0;0;0;0] 

% tol=0. 00005 
% gauss_seidel(A b,tol,x, N) 


Fig.(1.2) The GailSS-Seidel Matlab Program with its 
Command Window for RUNNING it. 
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function S0R(A, b,tol,x, N) 

%sor(A : b, tol : x ; N) solve iteratively a system of linear equations %whereby 
%A is the coefficient matrix, and b is the right-hand side column vector 
%tol : relative residual error tolerance for break condition 
%X : Nxl start vector (the initial guess) 

%N is the maximum number of iterations 
%The method implemented is that of Successive Over Relaxation 
%The starting vector is the null vector, but can be adjusted to one's needs 
%The iterative form is based on the SOR. transitioniteration matrix 


%Tw = inv(D-w*L)*((l-w)*D+w*U) and the constant vector 
% cw = w*inv(D-w*L)*b 


%The optimal parameter wis calculated using the spectral radius of the 
%Jacobi transition matrix Tj = inv(D)*(L+U) 

%The output is the solution vector X 


%splitting matrix A into the three matrices L, U and D 
D = diag(diag(A)); 

L = tril(-A,-1); 

U = triu(-A,l); 

Tj = inv(D)*(L+U); %Jacobi iteration matrix 
rho_Tj = max(abs(eig(Tj))); %spectral radius of Tj 
w= 1.25;%overrelaxation parameter 
dispCw - );disp(w); 

Tw= mv(D-w*L)*((l-w)*D+w*U); %SOR iteration matrix 
cw = w*inv(D-w*L)*b; %constant vector needed for iterations 


»A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1-1: 1-1 60 
0-2; 1 1 0 84 ; 0-1 -24700]; 

»N=32; 

»tol=0.00005 ; 

»x=[0:0;0;0:0]; 

»b=[U;3;4;5]; 

»SOR(A,b,tol,x,N) 
w= 1.2500 

The procedure was successful 
Condition ||x A (k+l)- x A (k)|| <tol 
was met after k iterations 15 


k= 1; 

while k<=N 

x(:Jc+l) = Tw*x(: Jc) + cw; 

B=norm(x(:,k+l)-x(:Jc)); 
if B<tol 

dispfThe procedure was successful') 

disp(Condition ||x A (k+l)- x A (k)|| <tol was met after k iterations') 
disp(k); dispfx = ');disp(x(:,k+l)); 


7.8597 

0.4229 

-0.0736 

-0.5406 

0.0106 


break 
end 

k=k+l; 
end 

if B >tol || k >=N 

dispCMaximum number of iterations reached without satisfying condition:') 
dispd|x A (k+l) - x A (k)|| <tol'); disp(tol); 
dispCPlease, examine the sequence of iterates') 

dispCIn case you observe convergence, then increase the maximum number of iterations'); disp(x); 
end 

% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&& 


Fig.(1.3) The SOR Matlab Program with its 

Command Window for RUNNING it. 


% A=[0.2 0.1 1 1 0 : 0.1 4 -1 1 -1: 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 
% N=30; 

% b=[l;2;3;4;5]; % x=[0;0;0;0;0]; 

% tol=0.00005 ; % SOR(A, b,toi,x, N) 


Now , let T w = (D -f wL) 1 [Cl — w] D — wU] 
and C w = wCD -f wL) _1 b 
Then the SOR technique has the iterative form 

= T w X ,k_1 ^ + for a constant w> 0, 
k=l,2,..., (1.14) 

By Eqn. (1.14) it is obvious that the method of SOR 

is an iterative technique that solves the left hand side of this 
expression for X'^' when the previous values for X k-i ' on the 
right hand side are computed .Analytically ,this may be 
written as: 

= (D +- wL ) -1 

(wb - [wU -I- (w - DDlX^), 
k=l,2 ,3... (1.15) 

By taking advantage of the triangular form of 
(D+ wL) it can proved that (D T wL ) -1 = D _1 . 

Using this, the elements of X^ can be computed sequentially 
using forward substitution as in Eqn.(1.3) in section l.land 

Eqn. (1.10) in section 1.2 above , i.e. 

(Id % (k-D 

Xj = U - wJXj 

+jr.( h i - Sj:,[ 3-ij x}'~ L ~ Sjwajj Xj* -1 ) , 

i=l,2,'...,n, k=0,l,.... (1.16) 

The choice of the relaxation factor w is not necessarily easy 
and depends upon the properties of the coefficient matrix .For 
positive -definite matrices it can be proved that 0<w<2will 
lead to convergence, but we are generally interested in faster 
convergence rather than just convergence. 


G. The Matlab Program for SOR iteration method 

The Matlab Program for SOR Method with it it’s 
Command Window is shown in Fig. (1.3) 


II. The conjugate gradient method 

It is used to solve the system of linear equations 
Ax = b (2.1) 

for the vector x where the known n-by-n matrix A is 
symmetric (i.e. A T = A), positive definite (i.e. x T Ax > 0 for all 
non-zero vectors x in R n ), and real, and b is known as well. 
We denote the unique solution of this system by x*. 

A. The conjugate gradient method as a direct method 

We say that two non-zero vectors u and v are 
conjugate (with respect to A) if 

•Tl 

11 Av = 0 . ( 2 . 2 ) 

Since A is symmetric and positive definite, the left-hand side 
defines an inner product 

(UV) A := (AU.V) =: (If, A r V) = (U, AV) = U r AV (2.3) 

Two vectors are conjugate if they are orthogonal with respect 
to this inner product. Being conjugate is a symmetric relation: 
if u is conjugate to v, then v is conjugate to u. (Note: This 
notion of conjugate is not related to the notion of complex 
conjugate.) 

Suppose that { } is a sequence of n mutually conjugate 
directions. Then the{ F k } form a basis of R n , so we can expand 
the solution x^of Ax = b in this basis: 

n 

i =1 

and we see that 

n 

b = AX. = ^ DCj A 

i = l 

The coefficients are given by 

n 

Ffb^F^AX.^^Ki F$AF* = cc K |£AIk (2.6) 

] = 1 


(2.4) 

(2.5) 


(because Vi ^ K, Pi,P K are mutually conjugate) 

K _ p K b _ < p K> b > _ < p K . b> Q 
K p gA p K ( p K ' p k)j4 IIPKII^ ' ' } 

This result is perhaps most transparent by considering the 
inner product defined above. 

This gives the following method for solving the equation Ax = 
b: find a sequence of n conjugate directions, and then 
compute the coefficients 

B. The conjugate gradient method as an iterative method 

The direct algorithm was then modified to obtain an 
algorithm which only requires storage of the last two residue 
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vectors and the last search direction, and only one matrix 
vector multiplication. 

The algorithm is detailed below for solving Ax = b 
where A is a real, symmetric, positive-definite matrix. The 
input vector x 0 can be an approximate initial solution or 0. 

:= b - Ax : 

Pu := r 0 
K := 0 


repeat 

k ' P^Pk 

- V A r + l T fX K 

r K + 1 r K ~ lX it- ^ Pk 

if r K +i is sufficiently small then exit loop 


T 


Pk 

^ + 1 := ^k+1 + f^Pk 
k = k + 1 


end repeat 
The result is 


This is the most commonly used algorithm. The same formula 
for f> K is also used in the Fletcher-Reeves nonlinear 
conjugate gradient method. 


C. The code in Matlab for the Conjugate 

Iterative Algorithm in Fiq.(1.4) & Fiq.(1.5) 

Alternatively another Matlab Code for Conjugate 
Gradient Algorithm is in Fiq.(1.5) 


D. Convergence properties of the conjugate gradient 
method 

The conjugate gradient method can theoretically be viewed 
as a direct method, as it produces the exact solution after a 
finite number of iterations, which is not larger than the size of 
the matrix, in the absence of round-off error. However, the 
conjugate gradient method is unstable with respect to even 
small perturbations, e.g., most directions are not in practice 
conjugate, and the exact solution is never obtained. 
Fortunately, the conjugate gradient method can be used as an 
iterative method as it provides monotonically improving 
approximations x k to the exact solution, which may reach 

the required tolerance after a relatively small (compared to the 
problem size) number of iterations. The improvement is 
typically linear and its speed is determined by the condition 
number K [A j of the system matrix A: the larger is K(A'), the 

slower the improvement. 

If K(A) is large, preconditioning is used to 


E. The preconditioned conjugate gradient method 

Preconditioning is necessary to ensure fast convergence 
of the conjugate gradient method. The preconditioned 
conjugate gradient method takes the following form: 
r 0 := b — Axn 
z c , := 

P» := 

K := 0 




repeat 
T 


P^A P k 


+ l Pk 

^ir+l := fit — Pk 

if r^ +1 is sufficiently small then exit loop end if 
Zk+1 : = 

T 

G . z k+t r k+i 

Pk ' — T _ 

Pit +1 := z k + t + ft it Pk 

k :=k+l 


end repeat 

The result is %- K +i 

The above formulation is equivalent to applying the 
conjugate gradient method without preconditioning to the 
system 

E- 1 A(E- 1 0 t x= E -1 b 
where EE^ = M and x = E T x 


The preconditioned matrix M has to be symmetric 
positive-definite and fixed, i.e., cannot change from iteration 
to iteration. If any of these assumptions on the preconditioned 
is violated, the behavior of the preconditioned conjugate 
gradient method may become unpredictable. 

An example of a commonly used preconditioned is the 
incomplete Cholesky factorization. 


F. The code in Matlab for the Preconditioned Conjugate 
Gradient Algorithm in Fig(1.6) 


G. The flexible preconditioned conjugate gradient 
method 

In numerically challenging applications 
sophisticated preconditioners are used, which may 
lead to variable preconditioning, changing between 
iterations. Even if the preconditioned matrix is 
symmetric positive-definite on every iteration, the 
fact that it may change makes the arguments above invalid, 
and in practical tests leads to a significant slowdown of the 
convergence of the algorithm presented above. Using the 
Polak-Ribiere formula 
„ *3r+ ifr'jr+i - 


replace the original system 
Ax — b = 0 with Af -1 (Ax — b) = 0 so that 

K gets smaller than K GO , see below 


instead of the Fletcher-Reeves formula 


Pk- 


r K+l 

Z K + l r K 


may dramatically improve the convergence in this case. This 

version of the preconditioned conjugate 

gradient method can be called flexible, as it allows 
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function [x, niter, flag] = CONJUGATE_GRAD(A, b, s, tol, maxiter) 
%CG Conjugate Gradients method 
°/oInput parameters 

°/oA : Symmetric, positive definite NxN matrix 
°/ob : Right-hand side Nxl column vector 
°/os : Nxl start vector (the initial guess) 

°/otol : relative residual error tolerance for break condition 
°/omaxiter : Maximum number of iterations to perform 

%Output parameters 

%x : Nxl solution vector 

%niter : Number of iterations performed 

%flag : 1 if convergence criteria specified by TOL could 

%notbe fulfilled within (the specified maximum Successful) 

%number of iterations, 0 otherwise 
x= s; % Set xO to the start vector s 
r= b - A*s; % Compute first residuum 
P = r; 


rho = r'*r; 

niter = 0; % Init counter for number of iterations 
flag = 0; % Init break flag 

%Computenorm of right-hand side to take relative residuum as 
%break condition 
normb = norm(b); 

if normb<eps % if the norm is vers close to zero, take the 
%absolute residuum instead as break condition 
%norm(r) >tol since the relative 

warning(['norm(b) is very close to zero, taking absolute as a break condition']); 


normb = 1; 
end 

while (norm(r)/normb>toI) 
a = A*p; 

alpha = rho/(a'*p); 
x = x+ alpha'p; 
r = r alpha'a; 
rhonew = r'*r; 
p = r + rho new/rho * p; 
rho = rhonew; 
niter = niter + 1; 

% if max. number of iterations 

if (niter — maxiter) 

flag =1; % is reached, break 


% Test break condition 


» A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 08 4 ; 0 -1 -2 4 700]; 
» maxiter=6; 

» b=[l;2;3;4;5]; 

» s=[0;0;0;0;0]; 

» tol=0.00005; 

» [x, niter, flag] = CONJUGATE_GRADIENT(A, b, s, tol, maxiter) 

x = 7.8597 0.4229 -0.0736 -0.5406 0.0106 
niter = 5 
flag = 0 


break 

end 

end 


% &&&&& ES THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&& 
% A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 

% maxiter=6; 

% b=[l;2;3;4;5]; Fig.(1.4) The CON JUGATE GRADIENT Matlab Program with its 

% s=[0;0;0;0;0]; Command Window for RUNNING it & below it the output 

% tol=0.00005; 

% [x, niter, flag] = CONJUGATE_GRADIENT(A, b, s, tol maxiter) 


particular, it converges not slower than the locally optimal 
steepest descent method. 


III. Conclusion 

When We ran each Matlab Code program and testified 
with a numerical example We observe that 
the non-stationary algorithms converge faster than the 
stationary algorithms. 

The reader can look and compare at the output result 
programs in Figures denoted above ((these Figures are 
Fig.(l) , Fig(2) , ....,Fig(5) and Fig(6) )) to ensure that 
non-stationary algorithms especially 

the Preconditioned Conjugate Algorithm approximates the 
solution accurately to five decimal 

places with the number of iterations N equals 5 which less 
compared to the other Algorithms tackled in this 
paper. 

The paper also analysis the Convergence properties of the 
conjugate gradient method and shows there are two methods, 
the direct or the iterative conjugate. 

Further, better convergence behavior can be achieved 
when using Polak-Ribiere formula which 
converges locally optimal not slower than the locally 
optimal steepest descent methed. 

Furthermore, these Algorithms can are 
recommended for similar situations which are arise in many 
important settings such as finite differences, 
finite element methods for solving partial differential 
equations and structural and circuit analysis. 


function [iteratioimo , x] = ccnjgrad{A ; b ; x) 
r=b-A*x; 


p=r; 

rsold=r'*r; 


^A=[Q.2 0.1 1 1 0 ;0 1 4 -1 1 - 1 : 1 -1 600 -2; 1 1 0 8 4 : 0 -1 -2 
4 700 ]; 


for 1=1:10000000 

Ap=A*p; 

alpha=rsold/ (p" * Ap) ; 
x=x— alpha*p; 
r=r- alpha :+: Ap; 


»b=|l;2;3;4;5]; 

»x=[0;0;0;0;0]; 

» [iterationno , x] = conj.grad{A,b ; x) 
iterationno =6 


rmew=r'*r; 
if sqit(rsnew)<5e- 10 
break; 

end 

p=r+rsnew/ rsold * p ; 
rsold=rsnew; 


x= 7.8597 
0.4229 
-0.073d 
-0.5406 
0.0106 


end 

iterationno=i 


Figil .5) The CON JUGATEGRADIENT Mafb b Program with 
its Command Window for RUNNING it 
& below it the output result 


% & IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM & 


% A=[0.2 0.1 1 1 0 ;0.1 4-1 1 -1; 1 -1 60 0-2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 
% b=[l;2;3;4;5]; 

% x=[0;0;0;0;0]; 

% [iterationno , x] = conjgrad{A,b,x) 


function Preconditioned_Conjugate(Ab,c : toLN) 

[njn]=size(A); 
invc=inv(c); 
x=zeros(n,l); 
r=b; 

p=invc*b; 
y^nvc*r; 
err(l)=norm(r); 
k=0; 

while noim(r) >tol&&k<N 
k=k+l; 
z=A*p; 

alpha=(y’*r)(p'*z); 
x=x-alpha*p; 
me w=r-alpha *z; 
ynew=invc *mew; 
beta=(ynew'*mewr) (y'*r); 
p=ynew+beta*p; 
r=mew; 
v=ynew; 
eir(k)=noim(r); 
end 

dispf x 3 '); disp(x'); 
dispf en=0; disp(eir); 

% &&&&& IN THE COMMAND \WNDOW TYPE TO RUN THE PROGRAM 


% A=[0J2 0.1 1 1 0 ; 0.1 4-1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 
% N=5 ; 

% b=[l;2;3;4;5]; 

% c=eye(5,5); 

% tol=<) .00005 ; 

% Preconditioned_Conjugate(A,b : c,toLN) 


Fig.(1.6) The Preconditioned Conjugate Matlab Program with 
Command Window for RUNNING it & below it the output result 


» A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1 -1) 1 -160 0 -2» 1 1 0 8 4 $ 0 -1 -2 4 700] 

>>N=5 ; 

» b=[l;2;3;4;5]; 

» c=eye(5,5); 

» tol=0. 00005 ; 

»Preconditioned_Conjugate(A,b,c,tol,N) 

x= 7.8597 0.4229 -0.0736 -0.5406 0.0106 

err= 7.5271 5.5600 0.7239 0.5572 0.0000 


for variable preconditioning. The implementation of 
the flexible version requires storing an extra vector. 

For a fixed preconditioned, = 0 so both 

formulas for are equivalent in exact arithmetic, i.e., 

without the round-off error. 

The mathematical explanation of the better 
convergence behavior of the method with the Polak-Ribiere 
formula is that the method is locally optimal in this case, in 
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